{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This script aims to extract model sources in a clear and informative format.  The script first shows what all the kinetic and thermo sources are in a model. Then it goes through each reaction and species to show their source and what the assigned uncertainties are.  This can be used with any RMG-generated CHEMKIN file that is annotated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from rmgpy.tools.uncertainty import Uncertainty, ThermoParameterUncertainty, KineticParameterUncertainty\n",
    "from IPython.display import display\n",
    "import copy\n",
    "import numpy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "chemFile = 'data/parseSource/chem_annotated.inp'\n",
    "dictFile = 'data/parseSource/species_dictionary.txt'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "uncertainty = Uncertainty(outputDirectory='testUncertainty')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "uncertainty.loadModel(chemFile, dictFile)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "uncertainty.loadDatabase()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "uncertainty.extractSourcesFromModel()\n",
    "uncertainty.compileAllSources()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print 'All Kinetic Sources'\n",
    "for sourceType in uncertainty.allKineticSources.keys():\n",
    "    if sourceType == 'Library':\n",
    "        print '============'\n",
    "        print 'Library kinetics'\n",
    "        print ''\n",
    "        print '\\tReactions: ', uncertainty.allKineticSources['Library']\n",
    "    elif sourceType == 'PDep':\n",
    "        print '============'\n",
    "        print 'PDep kinetics'\n",
    "        print ''\n",
    "        print '\\tReactions: ', uncertainty.allKineticSources['PDep']\n",
    "    elif sourceType == 'Rate Rules':\n",
    "        print '============'\n",
    "        print 'Rate rule kinetics'\n",
    "        print ''\n",
    "        for familyLabel, entries in uncertainty.allKineticSources['Rate Rules'].iteritems():\n",
    "            print '\\t', familyLabel\n",
    "            for entry in entries:\n",
    "                print '\\t\\t', entry\n",
    "    elif sourceType == 'Training':\n",
    "        print '============'\n",
    "        print 'Training reaction kinetics'\n",
    "        print ''\n",
    "        for familyLabel, entries in uncertainty.allKineticSources['Training'].iteritems():\n",
    "            print '\\t', familyLabel\n",
    "            for entry in entries:\n",
    "                print '\\t\\t', entry\n",
    "    else:\n",
    "        print sourceType\n",
    "        raise Exception('Kinetics source must be Library, PDep, Rate Rules, or Training')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print 'All Thermo Sources'\n",
    "for sourceType in uncertainty.allThermoSources.keys():\n",
    "    if sourceType == 'Library':\n",
    "        print '============'\n",
    "        print 'Library thermo'\n",
    "        print ''\n",
    "        print '\\tSpecies: ', uncertainty.allThermoSources['Library']\n",
    "    elif sourceType == 'QM':\n",
    "        print '============'\n",
    "        print 'QM thermo'\n",
    "        print ''\n",
    "        print '\\tSpecies: ', uncertainty.allThermoSources['QM']\n",
    "    elif sourceType == 'GAV':\n",
    "        print '============'\n",
    "        print 'Group additivity thermo'\n",
    "        print ''\n",
    "        for groupType, entries in uncertainty.allThermoSources['GAV'].iteritems():\n",
    "            print '\\t', groupType\n",
    "            for entry in entries:\n",
    "                print '\\t\\t', entry\n",
    "    else:\n",
    "        raise Exception('Thermo source must be GAV, QM, or Library')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Assign all the uncertainties using the Uncertainty() class function\n",
    "# ThermoParameterUncertainty and KineticParameterUncertainty classes may be customized and passed into this function\n",
    "# if non-default constants for constructing the uncertainties are desired\n",
    "uncertainty.assignParameterUncertainties()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "T = 623 # temperature in Kelvin for which to evaluate kinetics\n",
    "P = 1e5  # Pa "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "gParamEngine = ThermoParameterUncertainty()\n",
    "kParamEngine = KineticParameterUncertainty()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for rxn, source in uncertainty.reactionSourcesDict.iteritems():\n",
    "    print '======'\n",
    "    print rxn\n",
    "    display(rxn)\n",
    "    if 'Library' in source:\n",
    "        print 'Library reaction'\n",
    "        print source['Library']\n",
    "    elif 'PDep' in source:\n",
    "        print 'PDep reaction'\n",
    "        print source['PDep']\n",
    "    elif 'Rate Rules' in source:\n",
    "        print 'Rate rule estimate'\n",
    "        family = source['Rate Rules'][0]\n",
    "        sourceDict = source['Rate Rules'][1]\n",
    "        originalTemplate = sourceDict['template']\n",
    "        print '\\tFamily = ', family\n",
    "        print '\\tOriginal Template = ', [group.label for group in originalTemplate]\n",
    "        print '\\tExact = ', sourceDict['exact']\n",
    "        rules = sourceDict['rules']\n",
    "        training = sourceDict['training']\n",
    "        if rules:\n",
    "            print '\\tRate rule sources:'\n",
    "            for ruleEntry, weight in rules:\n",
    "                print '\\t\\t', ruleEntry, '=', weight\n",
    "        if training:\n",
    "            print '\\tTraining sources:'\n",
    "            for ruleEntry, trainingEntry, weight in training:\n",
    "                print '\\t\\t', ruleEntry , 'mapped to', trainingEntry , '=', weight\n",
    "    elif 'Training' in source:\n",
    "        print 'Training reaction'\n",
    "        family = source['Training'][0]\n",
    "        training = source['Training'][1]\n",
    "        print '\\t Family = ', family\n",
    "        print '\\t\\t', training\n",
    "\n",
    "\n",
    "\n",
    "    print ''\n",
    "    print 'Rate coefficient at {} K = {:.2e}'.format(T, rxn.kinetics.getRateCoefficient(T,P))\n",
    "\n",
    "    # Uncomment the following lines if you want to verify that the parsing has been performed correctly by\n",
    "    # checking the values for both the original and reconstructed kinetics\n",
    "\n",
    "#     print '---------'\n",
    "#     print 'Original kinetics:'\n",
    "#     print rxn.kinetics\n",
    "#     print ''\n",
    "#     print 'Reconstructed kinetics from parsing:'\n",
    "#     reconstructedKinetics=uncertainty.database.kinetics.reconstructKineticsFromSource(rxn,source,fixBarrierHeight=True)\n",
    "#     print reconstructedKinetics\n",
    "\n",
    "    rxnIndex = uncertainty.reactionList.index(rxn)\n",
    "    print 'Uncertainty dln(k) = ', uncertainty.kineticInputUncertainties[rxnIndex]\n",
    "    \n",
    "#     # Test that the partial uncertainty calculation is working\n",
    "#     dlnk = 0.0\n",
    "#     if 'Rate Rules' in source:\n",
    "#         family = source['Rate Rules'][0]\n",
    "#         sourceDict = source['Rate Rules'][1]\n",
    "#         rules = sourceDict['rules']\n",
    "#         training = sourceDict['training']\n",
    "#         for ruleEntry, weight in rules:\n",
    "#             dlnk += kParamEngine.getPartialUncertaintyValue(source, 'Rate Rules', corrParam=ruleEntry, corrFamily=family)\n",
    "#         for ruleEntry, trainingEntry, weight in training:\n",
    "#             dlnk += kParamEngine.getPartialUncertaintyValue(source, 'Rate Rules', corrParam=ruleEntry, corrFamily=family)\n",
    "#         dlnk += kParamEngine.getPartialUncertaintyValue(source, 'Estimation')\n",
    "#     elif 'PDep' in source:\n",
    "#         dlnk += kParamEngine.getPartialUncertaintyValue(source, 'PDep', source['PDep'])\n",
    "#     elif 'Library' in source:\n",
    "#         dlnk += kParamEngine.getPartialUncertaintyValue(source, 'Library', source['Library'])\n",
    "#     elif 'Training' in source:\n",
    "#         dlnk += kParamEngine.getPartialUncertaintyValue(source, 'Training', source['Training'])\n",
    "#     print 'Uncertainty dlnk calculated using sum of partial values = ', dlnk\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for species, source in uncertainty.speciesSourcesDict.iteritems():\n",
    "    print '=========='\n",
    "    print species\n",
    "    display(species)\n",
    "    if 'Library' in source:\n",
    "        print 'Thermo Library: ', source['Library']\n",
    "    if 'QM' in source:\n",
    "        print 'QM: ', source['QM']\n",
    "    if 'GAV' in source:\n",
    "        print 'Group additivity:'\n",
    "        for groupType, groupList in source['GAV'].iteritems():\n",
    "            print '\\t', groupType\n",
    "            for group, weight in groupList:\n",
    "                print '\\t\\t', group, '=', weight\n",
    "                \n",
    "    \n",
    "                \n",
    "    spcIndex = uncertainty.speciesList.index(species)    \n",
    "    print ''\n",
    "    print 'Uncertainty dG = ', uncertainty.thermoInputUncertainties[spcIndex], ' kcal/mol'\n",
    "    \n",
    "    \n",
    "    # Test that the partial uncertainty calculation is working\n",
    "    dG = 0.0\n",
    "    if 'Library' in source:\n",
    "        dG += gParamEngine.getPartialUncertaintyValue(source, 'Library', corrParam=source['Library'])\n",
    "    if 'QM' in source:\n",
    "        dG += gParamEngine.getPartialUncertaintyValue(source, 'QM',corrParam=source['QM'])\n",
    "    if 'GAV' in source:\n",
    "        for groupType, groupList in source['GAV'].iteritems():\n",
    "            for group, weight in groupList:\n",
    "                dG += gParamEngine.getPartialUncertaintyValue(source, 'GAV', group, groupType)\n",
    "        dG += gParamEngine.getPartialUncertaintyValue(source, 'Estimation')\n",
    "    print 'Uncertainty dG calculated using sum of partial values = ', dG, ' kcal/mol'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Assign correlated parameter uncertainties \n",
    "uncertainty.assignParameterUncertainties(correlated=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# See the thermo correlated parameter partial uncertainties\n",
    "uncertainty.thermoInputUncertainties"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# See the kinetics correlated parameter partial uncertainties\n",
    "uncertainty.kineticInputUncertainties"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
